/* Stellarium Web Engine - Copyright (c) 2022 - Stellarium Labs SRL
 *
 * This program is licensed under the terms of the GNU AGPL v3, or
 * alternatively under a commercial licence.
 *
 * The terms of the AGPL v3 license can be found in the main directory of this
 * repository.
 */


#include <stdio.h>
#include "erfa_wrap.h"
#include "swe.h"

// ut in MJD.
double deltat(double ut)
{
    double y, t;
    y = eraEpj(ERFA_DJM0, ut);

    // Stephenson Morrison Hohenker 2016 Delta T model
    static const double smh2016[54][6]={
//Row         Years                  Polynomial Coefficients
//  i      K_i    K_{i+1}        a_0         a_1         a_2         a_3
//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
/*  1 */  {-720.0,   400.0,   20550.593,  -21268.478,   11863.418,   -4541.129},
/*  2 */  { 400.0,  1000.0,    6604.404,   -5981.266,    -505.093,    1349.609},
/*  3 */  {1000.0,  1500.0,    1467.654,   -2452.187,    2460.927,   -1183.759},
/*  4 */  {1500.0,  1600.0,     292.635,    -216.322,     -43.614,      56.681},
/*  5 */  {1600.0,  1650.0,      89.380,     -66.754,      31.607,     -10.497},
/*  6 */  {1650.0,  1720.0,      43.736,     -49.043,       0.227,      15.811},
/*  7 */  {1720.0,  1800.0,      10.730,      -1.321,      62.250,     -52.946},
/*  8 */  {1800.0,  1810.0,      18.714,      -4.457,      -1.509,       2.507},
/*  9 */  {1810.0,  1820.0,      15.255,       0.046,       6.012,      -4.634},
/* 10 */  {1820.0,  1830.0,      16.679,      -1.831,      -7.889,       3.799},
/* 11 */  {1830.0,  1840.0,      10.758,      -6.211,       3.509,      -0.388},
/* 12 */  {1840.0,  1850.0,       7.668,      -0.357,       2.345,      -0.338},
/* 13 */  {1850.0,  1855.0,       9.317,       1.659,       0.332,      -0.932},
/* 14 */  {1855.0,  1860.0,      10.376,      -0.472,      -2.463,       1.596},
/* 15 */  {1860.0,  1865.0,       9.038,      -0.610,       2.325,      -2.497},
/* 16 */  {1865.0,  1870.0,       8.256,      -3.450,      -5.166,       2.729},
/* 17 */  {1870.0,  1875.0,       2.369,      -5.596,       3.020,      -0.919},
/* 18 */  {1875.0,  1880.0,      -1.126,      -2.312,       0.264,      -0.037},
/* 19 */  {1880.0,  1885.0,      -3.211,      -1.894,       0.154,       0.562},
/* 20 */  {1885.0,  1890.0,      -4.388,       0.101,       1.841,      -1.438},
/* 21 */  {1890.0,  1895.0,      -3.884,      -0.531,      -2.473,       1.870},
/* 22 */  {1895.0,  1900.0,      -5.017,       0.134,       3.138,      -0.232},
/* 23 */  {1900.0,  1905.0,      -1.977,       5.715,       2.443,      -1.257},
/* 24 */  {1905.0,  1910.0,       4.923,       6.828,      -1.329,       0.720},
/* 25 */  {1910.0,  1915.0,      11.142,       6.330,       0.831,      -0.825},
/* 26 */  {1915.0,  1920.0,      17.479,       5.518,      -1.643,       0.262},
/* 27 */  {1920.0,  1925.0,      21.617,       3.020,      -0.856,       0.008},
/* 28 */  {1925.0,  1930.0,      23.789,       1.333,      -0.831,       0.127},
/* 29 */  {1930.0,  1935.0,      24.418,       0.052,      -0.449,       0.142},
/* 30 */  {1935.0,  1940.0,      24.164,      -0.419,      -0.022,       0.702},
/* 31 */  {1940.0,  1945.0,      24.426,       1.645,       2.086,      -1.106},
/* 32 */  {1945.0,  1950.0,      27.050,       2.499,      -1.232,       0.614},
/* 33 */  {1950.0,  1953.0,      28.932,       1.127,       0.220,      -0.277},
/* 34 */  {1953.0,  1956.0,      30.002,       0.737,      -0.610,       0.631},
/* 35 */  {1956.0,  1959.0,      30.760,       1.409,       1.282,      -0.799},
/* 36 */  {1959.0,  1962.0,      32.652,       1.577,      -1.115,       0.507},
/* 37 */  {1962.0,  1965.0,      33.621,       0.868,       0.406,       0.199},
/* 38 */  {1965.0,  1968.0,      35.093,       2.275,       1.002,      -0.414},
/* 39 */  {1968.0,  1971.0,      37.956,       3.035,      -0.242,       0.202},
/* 40 */  {1971.0,  1974.0,      40.951,       3.157,       0.364,      -0.229},
/* 41 */  {1974.0,  1977.0,      44.244,       3.198,      -0.323,       0.172},
/* 42 */  {1977.0,  1980.0,      47.291,       3.069,       0.193,      -0.192},
/* 43 */  {1980.0,  1983.0,      50.361,       2.878,      -0.384,       0.081},
/* 44 */  {1983.0,  1986.0,      52.936,       2.354,      -0.140,      -0.166},
/* 45 */  {1986.0,  1989.0,      54.984,       1.577,      -0.637,       0.448},
/* 46 */  {1989.0,  1992.0,      56.373,       1.649,       0.709,      -0.277},
/* 47 */  {1992.0,  1995.0,      58.453,       2.235,      -0.122,       0.111},
/* 48 */  {1995.0,  1998.0,      60.677,       2.324,       0.212,      -0.315},
/* 49 */  {1998.0,  2001.0,      62.899,       1.804,      -0.732,       0.112},
/* 50 */  {2001.0,  2004.0,      64.082,       0.675,      -0.396,       0.193},
/* 51 */  {2004.0,  2007.0,      64.555,       0.463,       0.184,      -0.008},
/* 52 */  {2007.0,  2010.0,      65.194,       0.809,       0.161,      -0.101},
/* 53 */  {2010.0,  2013.0,      66.063,       0.828,      -0.142,       0.168},
/* 54 */  {2013.0,  2016.0,      66.917,       1.046,       0.360,      -0.282}
    };

    if (y<-720.)
    {
        double fact=(y-1825.0)/100.;
        return -320.0+32.5*fact*fact;
    }
    if (y>2016.)
    {
        // Fabien: use part of Espenak & Meeus (2006), shifted to avoid
        // discontinuity in 2016
        t = y - 2000;
        double dt2016 = 62.92 + 0.32217 * 16 + 0.005589 * 16 * 16;
        return 62.92 + 0.32217 * t + 0.005589 * t * t - (dt2016 - 68.1024);
    }
    int i=0;
    while (smh2016[i][1]<y) i++;
    assert(i<54);
    t=(y-smh2016[i][0]) / (smh2016[i][1]-smh2016[i][0]);
    return ((smh2016[i][5]*t + smh2016[i][4])*t
        + smh2016[i][3])*t + smh2016[i][2];
}

/******** TESTS ***********************************************************/

#if COMPILE_TESTS


static void test_deltat(void)
{
     // From United States Naval Observatory table
    typedef struct deltat_test {
        int y;
        int m;
        int d;
        double deltat;
    } deltat_test_t;
    deltat_test_t usno_data[] = {
        {1980,  1,  1,  50.5387},
        {2000,  1,  1,  63.8285},
        {2010,  1,  1,  66.0699},
        {2017,  1,  1,  68.5927},
        {2018,  1,  1,  68.9677},
    };

    for (int i = 0; i < ARRAY_SIZE(usno_data); ++i) {
        double djm0, djm;
        deltat_test_t* testd = &usno_data[i];
        eraCal2jd(testd->y, testd->m, testd->d, &djm0, &djm);
        double dt = deltat(djm);
        bool res =tests_compare_time(dt / ERFA_DAYSEC,
                                  testd->deltat / ERFA_DAYSEC, 200);
        if (!res) {
            LOG_E("Delta T mismatch for year %d", testd->y);
            assert(0);
        }
    }
}

TEST_REGISTER(NULL, test_deltat, TEST_AUTO)

#endif
